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Abstract 

In a hybrid plasma solver (particle ions, fluid mass-less electrons) re- 
gions of vacuum, or very low charge density, can cause problems since the 
evaluation of the electric field involves division by charge density. This 
causes large electric fields in low density regions that can lead to numer- 
ical instabilities. Here we propose a self consistent handling of vacuum 
regions for hybrid solvers. 

Vacuum regions can be considered having infinite resistivity, and in 
this limit Faraday's law approaches a magnetic diffusion equation. We 
describe an algorithm that solves such a diffusion equation in regions 
with charge density below a threshold value. We also present an imple- 
mentation of this algorithm in a hybrid plasma solver, and an application 
to the interaction between the Moon and the solar wind. 

We also discuss the implementation of hyperresistivity for smoothing 
the electric field in a PIC solver. 



1 The hybrid equations 

In the hybrid approximation, ions are treated as particles, and electrons as a 
masslcss fluid. In what follows we use SI units. The trajectory of an ion, r(t) 
and v(t), with charge q and mass m, is computed from the Lorentz force, 

dr dv q , 

— =v, - = -E + vxB, 
dt dt m 

where E = E(r, t) is the electric field, and B = B(r, t) is the magnetic field. 
From now on we do not write out the dependence on r and t. The electric field 
is given by 

E = — (-J/ x B + /i^ 1 (V x B) x B - Vp e ) + —V x B, (1) 
pi Po 

where pi is the ion charge density, J/ is the ion current density, p e is the electron 
pressure, r\ is the resistivity, and no = 4it ■ 10~ 7 is the magnetic constant. Then 
Faraday's law is used to advance the magnetic field in time, 

- = -VxE. (2) 
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Note that the unknowns are the position and velocity of the ions, and the 
magnetic field on a grid, not the electric field, since it can always be computed 
from Further details on the hybrid model used here, and the discretization, 
can be found in Holmstrom ( 2011a|b ). 



2 Vacuum regions 

In regions of low ion charge density, pi, the hybrid method can have numerical 
problems. We see from ([!]) that the electric field computation involves a division 
by pi. In what follows we will not write out the / subscript, i.e. p — pj. Thus, 
in low density regions we will have large electric fields, and in the limit of zero 
charge density, the electric field magnitude will tend to infinity. This can lead 
to numerical instabilities, due to large gradients in the electric field, and due to 
large accelerations of ions. The solution quickly becomes unstable. 

We can either have a region where p = physically, like inside a resistive 
obstacle, or due to the statistical nature of particle in cell solvers, there can be 
regions where we simply have no macro particles, or very few. From now on we 
denote all the different cases of low density regions as vacuum regions. 

Many ad hoc solutions have been proposed to handle vacuum regions. One 
way is to set a minimum allowed ion charge density, e.g., if p is below a threshold 
value in a cell, p is set to that value for the cell. A similar solution is to set a 
maximum allowed value for the electric field. One can also introduce new ion 
sources that were not part of the original problem, to keep p large enough in all 
cells. If we have an absorbing obstacle we can instead of removing absorbed ions 
reduce their weight gradually over time toward zero. None of these solutions 
solve the original problem. In the case of threshold values we do not get the 
solution to the hybrid equations, and in the case of artificial sources or losses, we 
are solving the equations self consistently, but we are solving a different physical 
problem. 

However, a self consistent, physically correct way of handling the problem 



of vacuum regions was proposed a long time ago. Hewett (19801 noted that 
vacuum regions can be viewed as having infinite resistivity, and an algorithm 
can be devised where the resistivity in is set to a large value in regions of low 



density. Harned ( 1982 ) also used a similar idea and solved a Laplace equation 
for the electric field in vacuum regions. However, solving Laplace equation 
over a complicated region that is also changing over time is not an easy task. 
Especially if one wants to do it on a parallel computer, since solving Laplace 
equation on a grid involves solving a system of linear equations. 

3 Solving a diffusion equation in vacuum regions 



Building on the idea of Hewett ( 1980[ ) of having a large resistivity in low density 



regions, let us see what a large resistivity implies. 

Let the resistivity be variable in space, and time, n — Tj(r,t). If we assume 
that the resistive term dominate in the expression (JlJ for the electric field, then 
Faraday's law ^ becomes 

f =-Vx(^VxB 
at \po 
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For a constant resistivity we have that 

"-V 2 B. (3) 



dt p 

This is a diffusion equation for the magnetic held, and the steady state solution 
will be a solution to the Laplace equation V 2 B = 0. 



4 Time step limits for the diffusion equation 

The diffusion equation for the magnetic field ([3]) is similar to the heat equation 
and gives a time step limit for stability of 

At<^, 

where At is the time step for an explicit time integrator. For large resistivities 
this will set the limit for the allowed length of the time step, e.g., if we increase 
rj by 10, At needs to be decreased by 10 for stability. 



5 Algorithm 

The details of the algorithm for handling low density regions is as follows. We 
assume that we have a resistivity, rj = rj(r,t), that is defined in all regions of 
the solution domain. A cell is denoted a vacuum cell if p < p v for the cell. The 
electric field in a vacuum cell is computed by ([I]), with 1/p = 0, i.e. only the 
resistive term is non-zero. If the vacuum cell is outside any obstacle, then we 
set the resistivity to a vacuum resistivity, rj = rj v . If the vacuum cell is in an 
obstacle, we keep the original resistivity. 

The advantage of such an approach is that vacuum regions and regions of 
a resistive object can be handled by the original solver. The time advance of 
Faraday's law ([2| is the same for all cells. We only compute the electric field 
differently in different cells. 



6 Hyperresistivity 



The leapfrog time stepping scheme that we use to discretize Faraday's law in 
time has no numerical diffusivity. To stabilize the computations one may have 
to smooth the solution. Here we choose to use a hyperresistivity term in the 
electric field computation. 



Hyperresistivity of different orders were introduced by Maron ( 2008 ) for an 
MHD solver. Here we introduce a hyperresistivity, T)h, and the expression for 
the electric field (JlJ becomes 



E 



Pi 



J/ X B + J X B — Vp e ) + T]J - %V 2 J, 



(4) 



where the current, J^^^VxB. The advantage of hyperresistivity compared 
to smoothing, that is often implemented in hybrid solvers, is that the higher or- 
der derivatives of the Laplacian diffuse high frequency structures, while preserv- 



ing lower frequency ones (Maron 2008). Also, the addition of hyperresistivity 



3 



can often increase the maximum stable time step. In our implementation, the 
Laplacian 

'd 2 j x d 2 j v d 2 j z 



2 T _ / " "Jx <-> <Jy 



V 3 ' dx 2 dy 2 dz l 

is discretized using standard second order finite difference stencils. The hyper- 
resistive term is present in all parts of the domain, also in the vacuum regions, 
and in the examples that follows, we have used a value of r\h — 5 • 10 14 . 

7 An application: The Lunar plasma wake 

A good example where vacuum and low density regions occur is the plasma wake 
behind the Moon, formed by the absorption on the dayside of the impinging solar 
wind plasma. More details on a hybrid model of the interaction between the 



Moon and the solar wind can be found in Holmstrom (2012). The simulation 
parameters used here are similar to those in that work, with an interplanetary 
magnetic field (IMF) at a 45° angle to the solar wind flow. 

The different resistivity regions at t — 30 s are shown in Fig. [T] Here rj is 
defined to be 0, except for the interior of the Moon (a sphere of radius 1730 km) 
where r\ = 10 7 57m. Cells with a relative ion charge density less than 0.0001, that 
are outside the obstacle (the Moon), was considered vacuum cell (p v = 0.0001), 
and there the resistivity was set to a vacuum resistivity r\ v = 10 s . We clearly 
see that the vacuum region is irregular. It also changes with time. 

By following the above approach to handle vacuum regions, we have intro- 
duced two new numerical parameters. A threshold ion charge density, p v , below 
which we consider a cell to be a vacuum cell; and a vacuum resistivity, r) v , that 
we set the resistivity to in such cells. Ideally, rj v should be as large as possible 
and p v should be as small as possible, for the solution to approach the solution 
to the original hybrid equations. What limits the value of r\ v is the time step 
limit discussed in Section |4j The computational time will increase in proportion 
to rj v for an explicit time integrator. The effect of different minimum density, 
p v , are shown in Fig. [2j We see that p v has to be small enough to get a smooth 
transition between vacuum and non-vacuum regions. The solution converges 
as p v is decreased, as can be seen especially in the central wake region. We 
can note that the relative p v is as small as one part in 10000 for the converged 
solution. This can be contrasted with the approach of just setting a minimum 
charge density. Then one typically has to choose a value of the relative den- 
sity on the order of a few %, otherwise numerical instabilities will develop. In 
contrast, the magnetic diffusion equation has to be applied only in cells of very 
low density. The method is also stable in the presence of large gradients in the 
resistivity. The example shown in Fig. [T] and [2] has a step change in resistivity, 
from to 10 7 at the Lunar surface. 



8 Summary and conclusions 

We have shown that vacuum regions in a plasma particle in cell solver can be 
handled by setting a high resistivity in those regions. This leads to the solution 
of a magnetic diffusion equation in such regions. A minimum charge density 
parameter decide what cells are vacuum, where the magnetic diffusion equation 
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Figure 1: Resistivity in different regions of the simulation domain at t — 30 s. A cut 
in the plane of the IMF. The solar wind flows in from the left. The blue region is the 
interior of the Moon with r) = 10 7 firm Red shows solar wind plasma with r\ — 0, and 
in yellow is the wake region with r/ v = 10 8 , where p < p v — 0.0001. 




Figure 2: The effect of different minimum charge densities, p v . Shown is the com- 
ponent of the magnetic field that is perpendicular to the plane containing the IMF. 
The minimum charge density, p v , decrease from top to bottom, with values of 0.01, 
0.001, and 0.0001. The cuts are planes that contain the IMF, and the black lines are 
the magnetic field lines. The solar wind flows in from the left. 
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should be advanced in time. It is also possible to include arbitrary resistive 
obstacles. The algorithm was exemplified by modeling the plasma interaction 
between the solar wind and the Moon. 

The advantage of the method is that it handles vacuum regions and resistive 
obstacles in a self consistent manner. Also, the magnetic diffusion equation 
only has to be applied in very low density cells, e.g., a relative density of 0.0001, 
compared to a reference density. The method seems stable to discontinuities in 
the resistivity. 

A disadvantage is that the time advance of the magnetic diffusion equation 
requires a small time step. This is however needed anyway for some problems, 
e.g., if we model the Moon and want to include crustal magnetic anomalies. A 
solution for future study would be to use an implicit time integrator to advance 
Faraday's law. 
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